
rm(list=ls())
load("3_MSAR_PA.Rdata")
library(coda)
library(ggplot2)
library(foreign)
library(grid)
class(jags.parsamples) <- "mcmc.list"

# Replicates MCMC figures appendix
pdf("AppendixB_coda1.pdf",height=6,width=8)
plot(jags.parsamples[,'sigma.a'],main="Sigma Epsilon")
dev.off()
pdf("AppendixB_coda2.pdf",height=6,width=8)
plot(jags.parsamples[,'sigma.b'],main="Sigma Nu")
dev.off()
pdf("AppendixB_coda3.pdf",height=6,width=8)
plot(jags.parsamples[,'phi[1]'],main="Beta Macro Lag")
dev.off()
pdf("AppendixB_coda4.pdf",height=6,width=8)
plot(jags.parsamples[,'phi[3]'],main="Beta ICS")
dev.off()
pdf("AppendixB_coda5.pdf",height=6,width=8)
plot(jags.parsamples[,'phi[4]'],main="Beta Pres Approval")
dev.off()
pdf("AppendixB_coda6.pdf",height=6,width=8)
plot(jags.parsamples[,'epsilon'],main="Pi")
dev.off()
pdf("AppendixB_coda7.pdf",height=6,width=8)
plot(jags.parsamples[,'mu[120]'],main="Alpha_120")
dev.off()

summary(statadata)
nrow(statadata)


# Remove 2011-12. for future change
#statadata <- statadata[1:(nrow(statadata)-8),]


epsilon1 <- jags.parsamples[,241][[1]]
epsilon2 <- jags.parsamples[,241][[2]]
epsilon3 <- jags.parsamples[,241][[3]]
combepsilon <- c(epsilon1,epsilon2,epsilon3)
combphi1 <- c(jags.parsamples[,482][[1]],jags.parsamples[,482][[2]],jags.parsamples[,482][[3]])
pdf("AppendixB_phieps.pdf",height=6,width=8)
scatter.smooth(combepsilon,combphi1, lpars=list(col="red",lwd=2),ylab="AR1 Coefficient Samples",xlab="Pi Samples")
dev.off()


phi <- rep(NA,4)
for (i in 482:485){
     phi[(i-481)] <- mean(unlist(lapply(jags.parsamples[,i],mean)))
}

delta <- rep(NA,240)
mu <- rep(NA,240)
ypred <- rep(NA,240)
ytilde <- rep(NA,240)

delta <- array(NA,dim=c(240,12000))
mu <- array(NA,dim=c(240,12000))
ypred <- array(NA,dim=c(240,12000))
ytilde <- array(NA,dim=c(240,12000))
for (i in 1:240){
    delta[i,] <- unlist(jags.parsamples[,i])
    mu[i,] <- unlist(jags.parsamples[,(240+1+i)])
    ypred[i,] <- unlist(jags.parsamples[,(487+i)])
    ytilde[i,] <- (mu[i,] + unlist(jags.parsamples[,483])*mean(statadata$party) + unlist(jags.parsamples[,484])*mean(statadata$partyics) + unlist(jags.parsamples[,485])*mean(statadata$partyapprove))/(1-unlist(jags.parsamples[,482]))
}

meandelta <- apply(delta,c(1),mean)
meanytilde <- apply(ytilde,c(1),mean)
meanmu <- apply(mu,c(1),mean)
meanypred <- apply(ypred,c(1),mean)

meanytilde2 <- (meanmu + phi[2]*mean(statadata$party) + phi[3]*mean(statadata$partyics) + phi[4]*mean(statadata$partyapprove))/(1-phi[1])

mean(meanytilde)
mean(meanytilde2)
mean(statadata$adjusted_gallup)

plot(meanytilde,type="l")
plot(statadata$adjusted_gallup)
lines(meanytilde)
lines(meanytilde2)
lines(meanypred,lty=2)
plot(statadata$date,meandelta)

sum((meanytilde-mean(statadata$adjusted_gallup))^2)/sum((statadata$adjusted_gallup-mean(statadata$adjusted_gallup))^2)
sum((meanytilde2-mean(statadata$adjusted_gallup))^2)/sum((statadata$adjusted_gallup-mean(statadata$adjusted_gallup))^2)
sum((meanypred-mean(statadata$adjusted_gallup))^2)/sum((statadata$adjusted_gallup-mean(statadata$adjusted_gallup))^2)
res <- statadata$adjusted_gallup - meanypred
sqrt(mean(res^2))

#summary(ols.fit)
#sqrt(mean((statadata$adjusted_gallup-ols.fit$fit)^2))


statadata$mu <- meanmu
statadata$ytilde <- meanytilde2
statadata$delta <- meandelta
statadata$ypred <- meanypred

mudiff <- mu[2:240,]-mu[1:239,]
mudiff <- mudiff*100
sqmudiff <- mudiff^2
absmudiff <- abs(mudiff)

futurevar <- array(NA,dim=c(240,12000))
futureabsdev <- array(NA,dim=c(240,12000))

for (i in 1:238){
    if (i<233){
        futureabsdev[i+1,] <- apply(absmudiff[i:(i+7),],c(2),sum)
        futurevar[i+1,] <- apply(sqmudiff[i:(i+7),],c(2),sum)/8
    }    else {
        futureabsdev[i+1,] <- apply(absmudiff[i:239,],c(2),sum)
        futurevar[i+1,] <- apply(sqmudiff[i:239,],c(2),sum)/(240-i)
    }
}

futurevar[1,] <- apply(sqmudiff[1:6,],c(2),sum)/7
futureabsdev[1,] <- apply(absmudiff[1:6,],c(2),sum)
futurevar[240,] <- sqmudiff[239,]
futureabsdev[240,] <- absmudiff[239,]

avgfvar <- apply(futurevar,c(1),mean)
avgabsdev <- apply(futureabsdev,c(1),mean)

avgfvartau <-  1/apply(futurevar,c(1),var,na.rm=TRUE)
lbfvar <- apply(futurevar,c(1),quantile,probs=c(.05),na.rm=TRUE)
ubfvar <- apply(futurevar,c(1),quantile,probs=c(.95),na.rm=TRUE)

futuresd <- sqrt(futurevar)
avgfsd <-  apply(futuresd,c(1),mean)


avgfsdtau <-  1/apply(futuresd,c(1),var,na.rm=TRUE)
lb5fsd <- apply(futuresd,c(1),quantile,probs=c(.05),na.rm=TRUE)
lb10fsd <- apply(futuresd,c(1),quantile,probs=c(.10),na.rm=TRUE)
ubfsd <- apply(futuresd,c(1),quantile,probs=c(.95),na.rm=TRUE)

plot(avgfsd,type="l",ylim=c(-0.5,2))
lines(lb5fsd,lty=3)
lines(lb10fsd,lty=3)

statadata$avgfvar <- avgfvar
statadata$avgfvartau <- avgfvartau
statadata$lbfvar <- lbfvar
statadata$avgfsd <- avgfsd
statadata$avgfsdtau <- avgfsdtau
statadata$lb5fsd <- lb5fsd
statadata$lb10fsd <- lb10fsd


combdelta <- rbind(delta,ifelse(runif(12000)<=combepsilon,1,0),ifelse(runif(12000)<=combepsilon,1,0),ifelse(runif(12000)<=combepsilon,1,0),ifelse(runif(12000)<=combepsilon,1,0),ifelse(runif(12000)<=combepsilon,1,0),ifelse(runif(12000)<=combepsilon,1,0),ifelse(runif(12000)<=combepsilon,1,0))

futureshift <- array(NA,dim=c(240,12000))

for (i in 1:240){
    futureshift[i,] <- ifelse(apply(combdelta[i:(i+7),],c(2),sum)>0,1,0)
}

avgfutureshift <- apply(futureshift,c(1),mean)
statadata$avgfutureshift <- avgfutureshift

elitepartychange <- read.dta("1_elitepartychange.dta")

statadata$fpartyideochange <- rep(elitepartychange$fpartyideochange[-(1:20)],each=4)
statadata$fpartyagendachange <- rep(elitepartychange$fpartyagendachange[-(1:20)],each=4)

write.dta(statadata,"4_conditionalupdatingmeasures.dta")


rm(list=ls())
                   
